#'
#' @TODO 绘制诺莫图，决策曲线，DCA曲线
#' @title ## 制诺莫图，决策曲线，DCA曲线
#' @description
#' @param Data 输入数据，*data.frame* ，需要包含`time`、`status`与其他临床特征等列
#' @param Features 纳入分析的临床特征,`Data`中的某几列
#' @param od 结果输出路径
#' @param TimeBreak 时间断点，默认1、3、5年
#' @return `NULL`
#' @usage
#' NomCalibrationDCA(Data = data_2, Features = c('Age', 'Gender', 'ECOG'), TimeBreak = c(1, 3, 5), od = '/Pub/Users/wangyk/project/tmp/12315/')
#'
#' @author *WYK*
#'
NomCalibrationDCA <- function(Data = NULL, Features = NULL, TimeBreak = c(1, 3, 5),
                              od = "./", TimeType = "OS") {
    dir_now <- paste0(od, "/")

    if (!dir.exists(dir_now)) {
        dir.create(dir_now, recursive = T)
    }

    # nomogram/------------
    library(rms)

    data_for_nom <- Data %>%
        select(all_of(Features), sample, time, status)

    dd <<- datadist(data_for_nom)
    options(datadist = "dd")
    options(na.action = "na.delete")
    # summary(data_for_nom$time)

    # Features <- c('Score', 'Age', 'Stage', 'margin_status')
    data_for_nom <- as.data.frame(data_for_nom)

    formula_used <- as.formula(str_c("Surv(time,status) ~", paste0(Features, collapse = "+")))

    coxm <- rms::cph(
        formula = formula_used, data = data_for_nom, dist = "lognormal",surv = T
        )



    surv_1 <- Survival(coxm) # 构建生存概率函数

    nom <- nomogram(coxm, fun = list(function(x) surv_1(365 * TimeBreak[1], x), function(x) {
        surv_1(365 *
            TimeBreak[2], x)
    }, function(x) surv_1(365 * TimeBreak[3], x)), funlabel = c(
        str_glue("{TimeBreak[1]}-year Survival"),
        str_glue("{TimeBreak[2]}-year Survival"), str_glue("{TimeBreak[3]}-year Survival")
    ))

    pdf(
        file = file.path(dir_now, "nomogram_classical.pdf"), width = 12, height = 7,
        pointsize = 12
    )
    # par(mar=c(1,1,1,1),cex=0.9)##mar 图形空白边界 cex 文本和符号大小
    plot(nom)
    dev.off()

    # calibration------- m大于四分之一样本数，小于三分之一样本数
    set.seed(1110)
    a <- nrow(data_for_nom) / 4 %>%
        as.integer()
    b <- nrow(data_for_nom) / 3 %>%
        as.integer()
    m_num <- sample(a:b, 1)

    f <- cph(
        formula = formula_used, data = data_for_nom, x = T, y = T, surv = T,
        time.inc = 365 * TimeBreak[1]
    )

    cal <- rms::calibrate(f,
        u = 365 * TimeBreak[1], cmethod = "KM", method = "crossvalidation",
        m = m_num, B = a
    )

    f1 <- cph(
        formula = formula_used, data = data_for_nom, x = T, y = T, surv = T,
        time.inc = 365 * TimeBreak[2]
    )
    cal1 <- rms::calibrate(f1,
        u = 365 * TimeBreak[2], cmethod = "KM", method = "crossvalidation",
        m = m_num, B = a
    )

    f2 <- cph(
        formula = formula_used, data = data_for_nom, x = T, y = T, surv = T,
        time.inc = 365 * TimeBreak[3]
    )
    cal2 <- rms::calibrate(f2,
        u = 365 * TimeBreak[3], cmethod = "KM", method = "crossvalidation",
        m = m_num, B = a
    )

    pdf(file = file.path(dir_now, "calibrate.pdf"), width = 5, height = 5, pointsize = 13)
    par(mgp = c(2.5, 1, 0))
    par(mar = c(4, 4, 1, 1))
    plot(cal,
        xlim = c(0, 1), ylim = c(0, 1), lwd = 2.5, lty = 1, errbar.col = c("#3B39D4"),
        cex.lab = 1.1, subtitles = F, riskdist = F, xlab = str_glue("Nomogram predicted {TimeType} Probability"),
        ylab = str_glue("Observed {TimeType} probability"), col = c("#3B39D4")
    )

    plot(cal1,
        lwd = 2.5, lty = 1, riskdist = F, errbar.col = c("#EA404D"), col = c("#EA404D"),
        add = T
    )

    plot(cal2,
        lwd = 2.5, lty = 1, riskdist = F, errbar.col = c("#411667"), col = c("#411667"),
        add = T
    )

    abline(0, 1, lty = 3, lwd = 2.5, col = c("#0D0D0D"))

    legend("bottomright", cex = 0.8, bg = NULL, legend = c(
        "Ideal", str_glue("{TimeBreak[1]}-year"),
        str_glue("{TimeBreak[2]}-year"), str_glue("{TimeBreak[3]}-year")
    ), col = c(
        "#0D0D0D",
        "#3B39D4", "#EA404D", "#411667"
    ), lwd = c(2, 2.5, 2.5, 2.5), lty = c(
        3, 1,
        1, 1
    ))
    dev.off()

    # save(list = 'data_for_nom', file =
    # str_glue('{dir_now}/data_for_nom.RData')) calibration------------

    # DCA决策曲线-----------
    # source('/home/innertech/RCodes/Project_wangyk/Codelib_YK/some/stdca.r')
    data_for_nom$years1 <- as.numeric(summary(survfit(f, newdata = data_for_nom),
        times = 365 * TimeBreak[1]
    )$surv)
    data_for_nom$years2 <- as.numeric(summary(survfit(f1, newdata = data_for_nom),
        times = 365 * TimeBreak[2]
    )$surv)
    data_for_nom$years3 <- as.numeric(summary(survfit(f2, newdata = data_for_nom),
        times = 365 * TimeBreak[3]
    )$surv)
    data_for_nom <- data.frame(data_for_nom, stringsAsFactors = F)

    km1 <- stdca(data = data_for_nom, outcome = "status", ttoutcome = "time", timepoint = 365 *
        TimeBreak[1], predictors = "years1", probability = T, xstop = 0.9)
    dev.off()
    km2 <- stdca(data = data_for_nom, outcome = "status", ttoutcome = "time", timepoint = 365 *
        TimeBreak[2], predictors = "years2", probability = T, xstop = 0.9)
    dev.off()
    km3 <- stdca(data = data_for_nom, outcome = "status", ttoutcome = "time", timepoint = 365 *
        TimeBreak[3], predictors = "years3", probability = T, xstop = 0.9)
    dev.off()

    result.path <- dir_now

    pdf(file = file.path(result.path, "Decision_curve_plot.pdf"), width = 6, height = 6)
    plot(km1$net.benefit$threshold, km1$net.benefit$all,
        type = "l", lwd = 1, lty = 2,
        col = "#3B39D4", las = 1, xlim = c(0, 0.6), ylim = c(-0.1, 0.5), xlab = "Threshold Probability",
        ylab = "Net Benefit"
    )
    lines(km1$net.benefit$threshold, km1$net.benefit$years1,
        type = "l", col = "#3B39D4",
        lwd = 2, lty = 1
    )

    lines(km2$net.benefit$threshold, km2$net.benefit$all,
        type = "l", col = "#EA404D",
        lwd = 1, lty = 2
    )
    lines(km2$net.benefit$threshold, km2$net.benefit$years2,
        type = "l", col = "#EA404D",
        lwd = 2, lty = 1
    )

    lines(km3$net.benefit$threshold, km3$net.benefit$all,
        type = "l", col = "#411667",
        lwd = 1, lty = 2
    )
    lines(km3$net.benefit$threshold, km3$net.benefit$years3,
        type = "l", col = "#411667",
        lwd = 2, lty = 1
    )
    abline(h = 0, col = "black", lwd = 2, lty = 1)

    legend(0.3, 0.5, cex = 0.8, legend = c(
        "none", str_glue("All ({TimeBreak[1]}-year {TimeType})"),
        str_glue("Nomogram ({TimeBreak[1]}-year {TimeType})"), str_glue("All ({TimeBreak[2]}-year {TimeType})"),
        str_glue("Nomogram ({TimeBreak[2]}-year {TimeType})"), str_glue("All ({TimeBreak[3]}-year {TimeType})"),
        str_glue("Nomogram ({TimeBreak[3]}-year {TimeType})")
    ), col = c(
        "black",
        "#3B39D4", "#3B39D4", "#EA404D", "#EA404D", "#411667", "#411667"
    ), text.col = c(
        "black",
        "#3B39D4", "#3B39D4", "#EA404D", "#EA404D", "#411667", "#411667"
    ), lwd = c(
        2,
        1, 2, 1, 2, 1, 2
    ), lty = c(1, 2, 1, 2, 1, 2, 1))
    dev.off()


    # # DCA curve---------- library(ggDCA) library(survival) library(rms)
    # library(patchwork)

    # formula_used <- as.formula(str_c('Surv(time,status) ~', paste0(Features,
    # collapse = '+'))) coxm <- cph(formula = formula_used, data =
    # data_for_nom, dist = 'lognormal', surv = T)

    # res <- map( c(Features), function(x) { formula_temp <-
    # as.formula(str_c('Surv(time,status) ~', x)) coxm <- cph(formula =
    # formula_temp, data = data_for_nom, dist = 'lognormal', surv = T) } )


    # Stage <- cph(Surv(time, status) ~ Stage, data = Data) Score <-
    # cph(Surv(time, status) ~ Score, data = Data) Nomogram <- cph(Surv(time,
    # status) ~ Stage + Score, data = Data)

    # labeltest1 <- dca(Stage, Score, Nomogram, model.names = c('Stage',
    # 'Score', 'Nomogram'), times = 365 * 1 ) labeltest2 <- dca(Stage, Score,
    # Nomogram, model.names = c('Stage', 'Score', 'Nomogram'), times = 365 * 2
    # ) labeltest3 <- dca(Stage, Score, Nomogram, model.names = c('Stage',
    # 'Score', 'Nomogram'), times = 365 * 3 )

    # # 作图 p1 <- ggplot(labeltest1, color = c('#41498e', '#E0AD16',
    # '#ED0000FF', '#42B540FF', '#26b9d3', 'black'), linetype = 1 ) +
    # labs(title = '1-year DCA') + theme( legend.position = 'top',
    # legend.justification = 'center', plot.title = element_text(face =
    # 'plain', size = 16, hjust = 0.5) ) plotout( p = p1, od = paste0(od,
    # '/DCA/'), name = '1year', w = 6, h = 6 )

    # p2 <- ggplot(labeltest2, color = c('#41498e', '#E0AD16', '#ED0000FF',
    # '#42B540FF', '#26b9d3', 'black'), linetype = 1 ) + labs(title = '2-year
    # DCA') + theme( legend.position = 'top', legend.justification = 'center',
    # plot.title = element_text(face = 'plain', size = 16, hjust = 0.5) )
    # plotout( p = p2, od = paste0(od, '/DCA/'), name = '2year', w = 6, h = 6 )

    # p3 <- ggplot(labeltest3, color = c('#41498e', '#E0AD16', '#ED0000FF',
    # '#42B540FF', '#26b9d3', 'black'), linetype = 1 ) + labs(title = '3-year
    # DCA') + theme( legend.position = 'top', legend.justification = 'center',
    # plot.title = element_text(face = 'plain', size = 16, hjust = 0.5) )
    # plotout( p = p3, od = paste0(od, '/DCA/'), name = '3year', w = 6, h = 6 )
}
# DCA决策曲线---------- stdca FUN
stdca <- function(data, outcome, ttoutcome, timepoint, predictors, xstart = 0.01,
                  xstop = 0.99, xby = 0.01, ymin = -0.05, probability = NULL, harm = NULL, graph = TRUE,
                  intervention = FALSE, interventionper = 100, smooth = FALSE, loess.span = 0.1,
                  cmprsk = FALSE) {

    # LOADING REQUIRED LIBRARIES
    require(survival)
    require(stats)

    # ONLY KEEPING COMPLETE CASES
    data <- data[complete.cases(data[c(outcome, ttoutcome, predictors)]), c(
        outcome,
        ttoutcome, predictors
    )]

    # outcome MUST BE CODED AS 0 AND 1
    if ((length(data[!(data[outcome] == 0 || data[outcome] == 1), outcome]) > 0) &
        cmprsk == FALSE) {
        stop("outcome must be coded as 0 and 1")
    }

    # data MUST BE A DATA FRAME
    if (class(data) != "data.frame") {
        stop("Input data must be class data.frame")
    }

    # xstart IS BETWEEN 0 AND 1
    if (xstart < 0 || xstart > 1) {
        stop("xstart must lie between 0 and 1")
    }

    # xstop IS BETWEEN 0 AND 1
    if (xstop < 0 || xstop > 1) {
        stop("xstop must lie between 0 and 1")
    }

    # xby IS BETWEEN 0 AND 1
    if (xby <= 0 || xby >= 1) {
        stop("xby must lie between 0 and 1")
    }

    # xstart IS BEFORE xstop
    if (xstart >= xstop) {
        stop("xstop must be larger than xstart")
    }

    # STORING THE NUMBER OF PREDICTORS SPECIFIED
    pred.n <- length(predictors)

    # IF probability SPECIFIED ENSURING THAT EACH PREDICTOR IS INDICATED AS A T
    # OR F
    if (length(probability) > 0 && pred.n != length(probability)) {
        stop("Number of probabilities specified must be the same as the number of predictors being checked.")
    }


    # IF harm SPECIFIED ENSURING THAT EACH PREDICTOR HAS A SPECIFIED HARM
    if (length(harm) > 0 && pred.n != length(harm)) {
        stop("Number of harms specified must be the same as the number of predictors being checked.")
    }

    # INITIALIZING DEFAULT VALUES FOR PROBABILITES AND HARMS IF NOT SPECIFIED
    if (length(harm) == 0) {
        harm <- rep(0, pred.n)
    }
    if (length(probability) == 0) {
        probability <- rep(TRUE, pred.n)
    }

    # THE PREDICTOR NAMES CANNOT BE EQUAL TO all OR none.
    if (length(predictors[predictors == "all" | predictors == "none"])) {
        stop("Prediction names cannot be equal to all or none.")
    }

    # CHECKING THAT EACH probability ELEMENT IS EQUAL TO T OR F, AND CHECKING
    # THAT PROBABILITIES ARE BETWEEN 0 and 1 IF NOT A PROB THEN CONVERTING WITH
    # A COX REGRESSION
    for (m in 1:pred.n) {
        if (probability[m] != TRUE & probability[m] != FALSE) {
            stop("Each element of probability vector must be TRUE or FALSE")
        }
        if (probability[m] == TRUE & (max(data[predictors[m]]) > 1 | min(data[predictors[m]]) <
            0)) {
            stop(paste(predictors[m], "must be between 0 and 1 OR sepcified as a non-probability in the probability option",
                sep = " "
            ))
        }
        if (probability[m] == FALSE) {
            model <- NULL
            pred <- NULL
            model <- coxph(Surv(data.matrix(data[ttoutcome]), data.matrix(data[outcome])) ~
                data.matrix(data[predictors[m]]))
            surv.data <- data.frame(0)
            pred <- data.frame(1 - c(summary(survfit(model, newdata = surv.data),
                time = timepoint
            )$surv))
            names(pred) <- predictors[m]
            data <- cbind(data[names(data) != predictors[m]], pred)
            print(paste(predictors[m], "converted to a probability with Cox regression. Due to linearity and proportional hazards assumption, miscalibration may occur.",
                sep = " "
            ))
        }
    }

    ######### CALCULATING NET BENEFIT #########
    N <- dim(data)[1]

    # getting the probability of the event for all subjects this is used for
    # the net benefit associated with treating all patients
    if (cmprsk == FALSE) {
        km.cuminc <- survfit(Surv(data.matrix(data[ttoutcome]), data.matrix(data[outcome])) ~
            1)
        pd <- 1 - summary(km.cuminc, times = timepoint)$surv
    } else {
        require(cmprsk)
        cr.cuminc <- cuminc(data[[ttoutcome]], data[[outcome]])
        pd <- timepoints(cr.cuminc, times = timepoint)$est[1]
    }

    # creating dataset that is one line per threshold for the treat all and
    # treat none strategies; CREATING DATAFRAME THAT IS ONE LINE PER THRESHOLD
    # PER all AND none STRATEGY
    nb <- data.frame(seq(from = xstart, to = xstop, by = xby))
    names(nb) <- "threshold"
    interv <- nb
    error <- NULL

    nb["all"] <- pd - (1 - pd) * nb$threshold / (1 - nb$threshold)
    nb["none"] <- 0

    # CYCLING THROUGH EACH PREDICTOR AND CALCULATING NET BENEFIT
    for (m in 1:pred.n) {
        nb[predictors[m]] <- NA

        for (t in 1:length(nb$threshold)) {
            # calculating number of true and false postives;
            px <- sum(data[predictors[m]] > nb$threshold[t]) / N

            if (px == 0) {
                error <- rbind(error, paste(predictors[m], ": No observations with risk greater than ",
                    nb$threshold[t] * 100, "%",
                    sep = ""
                ))
                break
            } else {
                # calculate risk using Kaplan Meier
                if (cmprsk == FALSE) {
                    km.cuminc <- survfit(Surv(data.matrix(data[data[predictors[m]] >
                        nb$threshold[t], ttoutcome]), data.matrix(data[data[predictors[m]] >
                        nb$threshold[t], outcome])) ~ 1)
                    pdgivenx <- (1 - summary(km.cuminc, times = timepoint)$surv)
                    if (length(pdgivenx) == 0) {
                        error <- rbind(error, paste(predictors[m], ": No observations with risk greater than ",
                            nb$threshold[t] * 100, "% that have followup through the timepoint selected",
                            sep = ""
                        ))
                        break
                    }
                    # calculate risk using competing risk
                } else {
                    cr.cuminc <- cuminc(
                        data[[ttoutcome]][data[[predictors[m]]] > nb$threshold[t]],
                        data[[outcome]][data[[predictors[m]]] > nb$threshold[t]]
                    )
                    pdgivenx <- timepoints(cr.cuminc, times = timepoint)$est[1]
                    if (is.na(pdgivenx)) {
                        error <- rbind(error, paste(predictors[m], ": No observations with risk greater than ",
                            nb$threshold[t] * 100, "% that have followup through the timepoint selected",
                            sep = ""
                        ))
                        break
                    }
                }
                # calculating NB based on calculated risk
                nb[t, predictors[m]] <- pdgivenx * px - (1 - pdgivenx) * px * nb$threshold[t] / (1 -
                    nb$threshold[t]) - harm[m]
            }
        }
        interv[predictors[m]] <- (nb[predictors[m]] - nb["all"]) * interventionper / (interv$threshold / (1 -
            interv$threshold))
    }
    if (length(error) > 0) {
        print(paste(error, ", and therefore net benefit not calculable in this range.",
            sep = ""
        ))
    }

    # CYCLING THROUGH EACH PREDICTOR AND SMOOTH NET BENEFIT AND INTERVENTIONS
    # AVOIDED
    for (m in 1:pred.n) {
        if (smooth == TRUE) {
            lws <- loess(data.matrix(nb[!is.na(nb[[predictors[m]]]), predictors[m]]) ~
                data.matrix(nb[!is.na(nb[[predictors[m]]]), "threshold"]), span = loess.span)
            nb[!is.na(nb[[predictors[m]]]), paste(predictors[m], "_sm", sep = "")] <- lws$fitted

            lws <- loess(data.matrix(interv[!is.na(nb[[predictors[m]]]), predictors[m]]) ~
                data.matrix(interv[!is.na(nb[[predictors[m]]]), "threshold"]), span = loess.span)
            interv[!is.na(nb[[predictors[m]]]), paste(predictors[m], "_sm", sep = "")] <- lws$fitted
        }
    }


    # PLOTTING GRAPH IF REQUESTED
    if (graph == TRUE) {
        require(graphics)

        # PLOTTING INTERVENTIONS AVOIDED IF REQUESTED
        if (intervention == TRUE) {
            # initialize the legend label, color, and width using the standard
            # specs of the none and all lines
            legendlabel <- NULL
            legendcolor <- NULL
            legendwidth <- NULL
            legendpattern <- NULL

            # getting maximum number of avoided interventions
            ymax <- max(interv[predictors], na.rm = TRUE)

            # INITIALIZING EMPTY PLOT WITH LABELS
            plot(
                x = nb$threshold, y = nb$all, type = "n", xlim = c(xstart, xstop),
                ylim = c(ymin, ymax), xlab = "Threshold probability", ylab = paste(
                    "Net reduction in interventions per",
                    interventionper, "patients"
                )
            )

            # PLOTTING INTERVENTIONS AVOIDED FOR EACH PREDICTOR
            for (m in 1:pred.n) {
                if (smooth == TRUE) {
                    lines(interv$threshold, data.matrix(interv[paste(predictors[m],
                        "_sm",
                        sep = ""
                    )]), col = m, lty = 2)
                } else {
                    lines(interv$threshold, data.matrix(interv[predictors[m]]),
                        col = m,
                        lty = 2
                    )
                }

                # adding each model to the legend
                legendlabel <- c(legendlabel, predictors[m])
                legendcolor <- c(legendcolor, m)
                legendwidth <- c(legendwidth, 1)
                legendpattern <- c(legendpattern, 2)
            }
        } else {
            # PLOTTING NET BENEFIT IF REQUESTED initialize the legend label,
            # color, and width using the standard specs of the none and all
            # lines
            legendlabel <- c("None", "All")
            legendcolor <- c(17, 8)
            legendwidth <- c(2, 2)
            legendpattern <- c(1, 1)

            # getting maximum net benefit
            ymax <- max(nb[names(nb) != "threshold"], na.rm = TRUE)

            # inializing new benfit plot with treat all option
            plot(x = nb$threshold, y = nb$all, type = "l", col = 8, lwd = 2, xlim = c(
                xstart,
                xstop
            ), ylim = c(ymin, ymax), xlab = "Threshold probability", ylab = "Net benefit")
            # adding treat none option
            lines(x = nb$threshold, y = nb$none, lwd = 2)
            # PLOTTING net benefit FOR EACH PREDICTOR
            for (m in 1:pred.n) {
                if (smooth == TRUE) {
                    lines(nb$threshold, data.matrix(nb[paste(predictors[m], "_sm",
                        sep = ""
                    )]), col = m, lty = 2)
                } else {
                    lines(nb$threshold, data.matrix(nb[predictors[m]]), col = m, lty = 2)
                }
                # adding each model to the legend
                legendlabel <- c(legendlabel, predictors[m])
                legendcolor <- c(legendcolor, m)
                legendwidth <- c(legendwidth, 1)
                legendpattern <- c(legendpattern, 2)
            }
        }
        # then add the legend
        legend("topright", legendlabel,
            cex = 0.8, col = legendcolor, lwd = legendwidth,
            lty = legendpattern
        )
    }

    # RETURNING RESULTS
    results <- list()
    results$N <- N
    results$predictors <- data.frame(cbind(predictors, harm, probability))
    names(results$predictors) <- c("predictor", "harm.applied", "probability")
    results$interventions.avoided.per <- interventionper
    results$net.benefit <- nb
    results$interventions.avoided <- interv
    return(results)
}

library(tidyverse)




# Data <- Data

# walk()

# Data %>% mutate_if(across(.cols = Features,.fns = as.factor))


# Data$A3_T <- as.numeric(factor(Data$A3_T)) Data$A4_N <-
# as.numeric(factor(Data$A4_N)) Data$A5_M <- as.numeric(factor(Data$A5_M))
# Data$A6_Stage <- as.numeric(factor(Data$A6_Stage))

# library(rmda) set.seed(1984) T.Stage.model <- decision_curve(OS~A3_T, data =
# Data, bootstraps = 500) set.seed(1984) N.Stage.model <-
# decision_curve(OS~A4_N, data = Data, bootstraps = 500)

# set.seed(1984) Riskscore.model <- decision_curve(OS~RiskScore, data = Data,
# bootstraps = 500) set.seed(1984) nomo.model <-
# decision_curve(OS~A3_T+A4_N+RiskScore, data = Data, bootstraps = 500)
# pdf('PDFs/tcga_cli_DCA.pdf', width = 5, height = 5) plot_decision_curve(
# list(T.Stage.model, N.Stage.model, Riskscore.model, nomo.model), curve.names
# = c('T.Stage model', 'N.Stage model', 'RiskScore model', 'Nomogram model'),
# xlim = c(0, 1), legend.position = 'topright', confidence.intervals = FALSE,
# col = color9) dev.off()
